///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file correlator.cpp
///
/// \brief   contains implementations of functions related to the correlator class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include <src/pcps.h>
#include <src/input.h>
#include <src/correlator.h>

////////////////////////////////////////////////////////////////////////////////////
/// \brief   Writes the correlator to a binary output stream
///
///   \param  ofs     the output filestream to write the correlator to
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
void pcps::Correlator<S>::dump_binary(std::ofstream & ofs) const {

  if ( _ti ) throw pcps::Exception("cannot write translationally invariant placeholder to file");
  if ( _site_vec.get() == 0 ) throw pcps::Exception("cannot write to binary a correlator with an empty _site_vec pointer");
  if (     _data.get() == 0 ) throw pcps::Exception("cannot write to binary a correlator with an empty _data pointer"    );
  const int nsites = _site_vec->size();
  ofs.write((const char *)&nsites, sizeof(int));
  ofs.write((const char *)&_fixed, sizeof(bool));
  ofs.write((const char *)&_configs_per_site, sizeof(int));
  ofs.write((const char *)&_nelements, sizeof(int));
  ofs.write((const char *)&_site_vec->at(0), nsites * sizeof(int));
  ofs.write((const char *)&_data->at(0), _nelements * sizeof(S));

}

// explicitly instantiate the template
template void pcps::Correlator< double               >::dump_binary(std::ofstream &) const;
template void pcps::Correlator< std::complex<double> >::dump_binary(std::ofstream &) const;

////////////////////////////////////////////////////////////////////////////////////
/// \brief   Reads the correlator from a binary input stream
///
///   \param  ifs     the input filestream to read the correlator from
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
void pcps::Correlator<S>::read_binary(std::ifstream & ifs) {

  int nsites;
  ifs.read((char *)&nsites, sizeof(int));
  ifs.read((char *)&_fixed, sizeof(bool));
  ifs.read((char *)&_configs_per_site, sizeof(int));
  ifs.read((char *)&_nelements, sizeof(int));
  _site_vec = boost::shared_ptr< std::vector<int> >( new std::vector<int>(nsites) );
  ifs.read((char *)&_site_vec->at(0), nsites * sizeof(int));
  _data = boost::shared_ptr< std::vector<S> >( new std::vector<S>(_nelements) );
  ifs.read((char *)&_data->at(0), _nelements * sizeof(S));

}

// explicitly instantiate the template
template void pcps::Correlator< double               >::read_binary(std::ifstream &);
template void pcps::Correlator< std::complex<double> >::read_binary(std::ifstream &);

////////////////////////////////////////////////////////////////////////////////////
/// \brief   Broadcasts the correlator to all processes in the given communicator.
///
///   \param  comm     mpi communicator used for the broadcast
///   \param  myrank   the mpi rank of this process
///   \param  root     the mpi rank of the process from which the correlator will be broadcast
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
void pcps::Correlator<S>::broadcast(const MPI::Comm & comm, const int myrank, const int root) {

  // be sure the correlator to be broadcast is initialized
  if (myrank == root) {
    if ( _ti ) throw pcps::Exception("cannot broadcast translationally invariant correlator");
    if ( _site_vec.get() == 0 ) throw pcps::Exception("cannot broadcast correlator with an empty _site_vec pointer");
    if (     _data.get() == 0 ) throw pcps::Exception("cannot broadcast correlator with an empty _data pointer"    );
  }

  int nsites;
  if (myrank == root)
    nsites = _site_vec->size();
  comm.Bcast(&nsites, 1, MPI::INT, root);
  comm.Bcast(&_fixed, sizeof(bool), MPI::CHAR, root);
  comm.Bcast(&_ti, sizeof(bool), MPI::CHAR, root);
  comm.Bcast(&_configs_per_site, 1, MPI::INT, root);
  comm.Bcast(&_nelements, 1, MPI::INT, root);
  if (myrank != root) {
    _ti_sites = boost::shared_ptr< std::vector<int> >();
    _site_vec = boost::shared_ptr< std::vector<int> >( new std::vector<int>(nsites) );
    _data = boost::shared_ptr< std::vector<S> >( new std::vector<S>(_nelements) );
  }
  if (nsites > 0)
    comm.Bcast(&_site_vec->at(0), nsites, MPI::INT, root);
  if (_nelements > 0)
    comm.Bcast(&_data->at(0), _nelements * sizeof(S), MPI::CHAR, root);

}

// explicitly instantiate the template
template void pcps::Correlator< double               >::broadcast(const MPI::Comm &, const int, const int);
template void pcps::Correlator< std::complex<double> >::broadcast(const MPI::Comm &, const int, const int);

////////////////////////////////////////////////////////////////////////////////////
/// \brief   Absorbs smaller correlators into larger correlators.
///
///   For each correlator whose sites are contained by other correlators, this
///   function absorbs that correlator's values into the containing correlators
///   and removes it from the vector of correlators.
///
///   \param  userinp      user input options
///   \param  correlators  vector of correlators that will be absorbed into one another
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
void pcps::absorb_correlators(const pcps::Input &userinp, std::vector< pcps::Correlator<S> > &correlators) {

  // sort the correlators ascending by size
  std::sort(correlators.begin(), correlators.end());

  // loop over the correlators, absorbing correlators when possible
  for (int j = 0; j < correlators.size(); ) {

    // get a reference to the smaller correlator
    pcps::Correlator<S> & small_corr = correlators.at(j);

    // initialize a vector of correlators that contain the small correlator
    std::vector< pcps::Correlator<S> > containing;

    for (int i = j+1; i < correlators.size(); i++) {

      // get a reference to the larger correlator
      pcps::Correlator<S> & big_corr = correlators.at(i);

      // do not absorb fixed or translationally invariant correlators
      if (big_corr.fixed() || small_corr.fixed() || big_corr.ti() || small_corr.ti()) continue;

      // determine whether the sites of small_corr are contained by the sites of big_corr
      bool contained = true;
      for (std::vector<int>::const_iterator site = small_corr.site_vec().begin(); site != small_corr.site_vec().end(); site++)
        if (std::count(big_corr.site_vec().begin(), big_corr.site_vec().end(), *site) == 0)
          contained = false;

      // if small_corr is contained by big_corr, add big_corr to the set of containing correlators
      if (contained)
        containing.push_back(big_corr);

    }

    // if any correlators contain small_corr, absorb it into them
    if (containing.size() > 0) {
      for (typename std::vector<pcps::Correlator<S> >::iterator it = containing.begin(); it != containing.end(); it++) {
        for (pcps::Cluster occs(userinp.nsites(), userinp.configs_per_site(), it->site_vec()); !occs.finished(); occs++)
          (*it)[occs] *= std::pow(small_corr[occs], pcps::unity<S>() / double(containing.size()));
      }
      correlators.erase( correlators.begin() + j );
    } else {
      j++;
    }

  }

}

// explicitly instantiate the template
template void pcps::absorb_correlators(const pcps::Input &, std::vector< pcps::Correlator< double               > > &);
template void pcps::absorb_correlators(const pcps::Input &, std::vector< pcps::Correlator< std::complex<double> > > &);

////////////////////////////////////////////////////////////////////////////////////
/// \brief   Create a correlator for a spin coherent state.
///
///   Constructs a spin coherent state with spin angle expectation values theta
///   and phi using a single correlator.
///
///   \param  site               the site on which to construct the spin coherent state
///   \param  configs_per_site   the number of configurations available to a single site
///   \param  theta              an angle, 0 <= theta < pi
///   \param  phi                an angle, 0 <= phi < 2*pi
///   \param  fixed              whether or not the correlator should be fixed
///
///   \return a one-site correlator representing the spin coherent state
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
pcps::Correlator<S> pcps::spin_coherent_correlator(const int site,
                                                        const int configs_per_site,
                                                        const double theta,
                                                        const double phi,
                                                        const bool fixed) {

  const double pi = 3.14159265358979;
  assert( theta >= 0 && theta <   pi );
  assert(   phi >= 0 && phi   < 2*pi );

  pcps::Correlator<S> retval(site, configs_per_site, pcps::unity<S>(), fixed);

  const int ss = configs_per_site - 1;

  const S ii = std::sqrt( -1.0 * pcps::unity<S>() );

  for (int p = 0; p <= ss; p++)
    retval[p] = std::sqrt(pcps::n_choose_m(ss,p)) * std::pow(std::cos(theta/2), ss-p) * std::pow(std::exp(ii * phi) * std::sin(theta/2), p);

  return retval;

}

// explicitly instantiate the template
template pcps::Correlator< double               > pcps::spin_coherent_correlator(const int, const int, const double, const double, const bool);
template pcps::Correlator< std::complex<double> > pcps::spin_coherent_correlator(const int, const int, const double, const double, const bool);

////////////////////////////////////////////////////////////////////////////////////
/// \brief   prints the data inside a correlator
///
///   \param  corr   the correlator to be printed
///
////////////////////////////////////////////////////////////////////////////////////

template <class S>
void pcps::print_correlator(const pcps::Correlator<S> & corr) {

  std::cout << "printing correlator on sites";
  for (std::vector<int>::const_iterator site = corr.site_vec().begin(); site != corr.site_vec().end(); site++)
    std::cout << boost::format(" %5i") % (*site);
  std::cout << ":" << std::endl;

  const int nsites = corr.site_vec().size();
  const long int nconfigs = pcps::pow((long int)(corr.configs_per_site()), nsites);

  std::vector<int> occs(nsites, 0);

  // loop over the correlator's sites' occupations
  for (long int config_index = 0; config_index < nconfigs; config_index++) {

    // compute compound index from the occupation pattern
    int cmpd_index = 0;
    for (int j = 0; j < nsites; j++)
      cmpd_index = corr.configs_per_site() * cmpd_index + occs[j];

    // print out the occupation pattern and the corresponding compound index and correlator value
    std::cout << "value for occupation";
    for (std::vector<int>::const_iterator occ = occs.begin(); occ != occs.end(); occ++)
      std::cout << boost::format(" %5i") % (*occ);
    std::cout << boost::format("   (cmpd = %9i)   is   %s    norm = %12.6f")
                 % cmpd_index % pcps::format_number<S>("%12.6f", corr[cmpd_index]) % std::abs(corr[cmpd_index]) << std::endl;

    // create the next occupation pattern
    int i = nsites - 1;
    occs[i] += 1;
    while (i > 0 && occs[i] == corr.configs_per_site()) {
      occs[i] = 0;
      i--;
      occs[i] += 1;
    }

  }

}

// explicitly instantiate the template
template void pcps::print_correlator(const pcps::Correlator< double               > &);
template void pcps::print_correlator(const pcps::Correlator< std::complex<double> > &);

